Redes Infecciosas

Mayela Fosado Mendoza

2023-01-24

Mayela Fosado

Introducción

Los modelos por compartimentos son una herramienta útil para modelar la transmisión de una enfermedad. Con la paquetería desolve en R, se pueden representar los modelos, añadiendo los parámetros de las enfermedades, de esta manera se puede obtener un plot que en el eje \(x\) muestre el tiempo y en el eje \(y\) la cantidad de personas que hay por estado (susceptible, infectado, recuperado, etc.) Sin embargo, a través del plot no se puede visualizar cómo es el contacto entre las personas y las conexiones que pueden existir entre ellas, es decir, que habría más personas con un mayor a riesgo a enfermedad si están en contacto con alguien infectado.

Este proyecto tiene el propósito de poder desarrollar una herramienta que permita visualizar ese contacto entre las personas, por medio de una función de redes. Se utilizaron dos modelos de compartimentos para realizarlo, el modelo SIR y el modelo SEIR. Ambos modelos fueron modelados por los dos métodos en este proyecto, y así mismo se realizaró una comparación entre ellos.

Librerías para utilizar:

library(igraph)
## 
## Adjuntando el paquete: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
library(ggplot2)
library(animation)
library(ranger)
library(viridis)
## Cargando paquete requerido: viridisLite
library(deSolve)

Ecuaciones

S: Susceptibles

E: Expuestos / Exposed

I: Infectados / Infected

R: Recuperados / Recovery

Modelo SIR

\[\dot S=-\beta SI\] \[\dot I=\beta SI - \gamma I\] \[\dot R=\gamma I\]

Modelo SEIR

\[\dot S=-\beta SI\] \[\dot E=\beta SI - \delta E\] \[\dot I=\delta E-\gamma I\] \[\dot R=\gamma I\]

Modelo por compartimentos

Modelo SIR

Modelo SEIR

Función de redes

Crear una red de origen

Adj <- barabasi.game(500, directed = FALSE)
## Warning: `barabasi.game()` was deprecated in igraph 2.0.0.
## ℹ Please use `sample_pa()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Función Infectious.Network para establecer los parámetros de la transmisión de una enfermedad, el nodo en el que comienza el infectado, los días que dura una persona en periodo infeccioso, etc.

Partes de la función x: Objeto listo para plotear de igraph (la red)

StateConfig: c(S, I, R) o c(S, E, I, R)

Colors: c(Color1, Color2, Color3) o c(Color1, Color2, Color3, Color4)

StateDays: ¿Cuántos días permanece una persona en ese estado? c(0, #, #, #)

StartNode: ¿Quién(es) va(n) a ser el primer infeccioso?

R0: Promedio de personas que van a ser contagiadas por una persona

Days: ¿Cuántos días quieres que itere? ###

length(StateConfig) == length(Colors) == length(StateDays)

Infectious.Network <- function(x, StateConfig, StateDays, StartNode, R0, Colors, Days) {  
  
  
  # Hacer una base de datos que contenga el conteo de los estados por dia
  if (length(StateConfig) == 3){
    Summary <- data.frame(0, 0, 0, 0)
    colnames(Summary) <- c("S", "I", "R", "Day")
  } else if (length(StateConfig) == 4) {
    Summary <- data.frame(0, 0, 0, 0, 0 )
    colnames(Summary) <- c("S", "E", "I", "R", "Day")
  }
  
  
  
  # Para los nombres de los plots
  if (length(StateConfig) == 3) {PlotName <- "SIR"
  Lol <- c("Susceptible", "Infected", "Recovered")
  } else if (length(StateConfig) == 4) {PlotName <- "SEIR"
  Lol <- c("Susceptible", "Exposed", "Infected", "Recovered")}
  d <- 1
  
  
  
  # Asignar a objetos la cantidad de dias que pasa una persona en un estado determinado
  if(length(StateDays) == 3) {S<-StateDays[1]
  I<-StateDays[2]
  R<-StateDays[3]} else if (length(StateDays) == 4){
    S<-StateDays[1]
    E<-StateDays[2]
    I<-StateDays[3]
    R<-StateDays[4]
  }
  
  
  
  # Para que quede fija la red
  LO <- layout_nicely(Adj)
  angle <- 7*100 * 20
  RotMat <- matrix(c(cos(angle),sin(angle),-sin(angle), cos(angle)), ncol=2)
  LO2 <- LO %*% RotMat
  
  
  
  
  Net <- x # Asignar a otro objeto
  AdjNet <- as_adj_list(Net) # Para obtener adyacencias
  
  
  # Añadir atributos a los vertices
  V(Net)$State <- c(rep("A", length(Net)))
  V(Net)$Days <- c(rep(0, length(Net)))  
  
  
  
  
  # Asociar colores a los estados en los plots
  Colors <- Colors
  my_color <- Colors[as.numeric(as.factor(V(Net)$State))]
  
  
  
  
  
  # Asignar estados a un orden alfabetico (porque igraph asigna los colores en orden alfabetico)
  for (i in 1:length(StateConfig)) {
    StateConfig[i] <- LETTERS[i]
  }
  
  
  
  ####################
  # HAZ UN PLOT AQUI #
  ####################
  my_color <- Colors[as.numeric(as.factor(V(Net)$State))]
  
  plot(Net, layout = LO2, vertex.color = my_color, vertex.size = degree(Adj, V(Adj), "in")*2, 
       vertex.label = "", edge.arrow.size = 0.2, edge.size = 19, 
       vertex.frame.color="white", edge.color = "#E5E3C9")
  
  legend(x=-1.5, y=-1.1, paste(Lol), pch=21, 
         col="lightgray", pt.bg = Colors, pt.cex=2, cex=1.5, bty="n", ncol=1, title = paste(PlotName))
  
  #######################################
  # Actualizar la base de datos por dia #
  #######################################
  if (length(StateConfig) == 3){
    new_row <- c(as.numeric(table(V(Net)$State == "A")["TRUE"]),
                 as.numeric(table(V(Net)$State == "B")["TRUE"]),
                 as.numeric(table(V(Net)$State == "C")["TRUE"]),
                 0)
    
    Summary <- rbind(Summary, new_row)
  } else if (length(StateConfig) == 4) {
    new_row <- c(as.numeric(table(V(Net)$State == "A")["TRUE"]),
                 as.numeric(table(V(Net)$State == "B")["TRUE"]),
                 as.numeric(table(V(Net)$State == "C")["TRUE"]),
                 as.numeric(table(V(Net)$State == "D")["TRUE"]),
                 0)
    
    Summary <- rbind(Summary, new_row)
  }
  
  
  
  # Asignar el primer infectado a la red
  V(Net)$Days[StartNode] <- 1
  if(length(StateConfig) == 3) {V(Net)$State[StartNode] <- "B"} else if (length(StateConfig) == 4) {V(Net)$State[StartNode] <- "C"}
  
  
  
  
  ####################
  # HAZ UN PLOT AQUI #
  ####################
  my_color <- Colors[as.numeric(as.factor(V(Net)$State))]
  
  plot(Net, layout = LO2, vertex.color = my_color, vertex.size = degree(Adj, V(Adj), "in")*2, 
       vertex.label = "", edge.arrow.size = 0.2, edge.size = 19, 
       vertex.frame.color="white", edge.color = "#E5E3C9")
  
  legend(x=-1.5, y=-1.1, paste(Lol), pch=21, 
         col="lightgray", pt.bg = Colors, pt.cex=2, cex=1.5, bty="n", ncol=1, title = paste(PlotName, "Day", 0))
  
  #######################################
  # Actualizar la base de datos por dia #
  #######################################
  if (length(StateConfig) == 3){
    new_row <- c(as.numeric(table(V(Net)$State == "A")["TRUE"]),
                 as.numeric(table(V(Net)$State == "B")["TRUE"]),
                 as.numeric(table(V(Net)$State == "C")["TRUE"]),
                 0)
    
    Summary <- rbind(Summary, new_row)
  } else if (length(StateConfig) == 4) {
    new_row <- c(as.numeric(table(V(Net)$State == "A")["TRUE"]),
                 as.numeric(table(V(Net)$State == "B")["TRUE"]),
                 as.numeric(table(V(Net)$State == "C")["TRUE"]),
                 as.numeric(table(V(Net)$State == "D")["TRUE"]),
                 0)
    
    Summary <- rbind(Summary, new_row)
  }
  
  
  
  for (i in 1:Days) {
    
    
    # Asignar los dias de los estados a nuevos objetos
    if(length(StateConfig) == 3) {S <- StateDays[1]
    I <- StateDays[2]
    R <- StateDays[3]
    } else if (length(StateConfig) == 4) {
      
      S <- StateDays[1]
      E <- StateDays[2]
      I <- StateDays[3]
      R <- StateDays[4]}
    
    
    
    
    
    # Transmision de la enfermedad
    if (length(StateConfig) == 3){ # primer if
      
      for (i in 1:length(V(Net)$State)) {
        if(V(Net)$State[i] == "B") {new_infected <- sample(as.numeric(AdjNet[[i]]), R0, replace = T)
        for (i in 1:length(new_infected)) {
          if (V(Net)$Days[new_infected[i]] == 0) {V(Net)$Days[new_infected[i]] <- 1}
        }
        }
      } 
      
      
    } else if (length(StateConfig) == 4) { # Segundo if 
      
      for (i in 1:length(V(Net)$State)) {
        if(V(Net)$State[i] == "C") {new_infected <- sample(as.numeric(AdjNet[[i]]), R0, replace = T)
        for (i in 1:length(new_infected)) {
          if (V(Net)$Days[new_infected[i]] == 0) {V(Net)$Days[new_infected[i]] <- 1}
        }
        }
      } 
      
      
    } # Aqui termina el segundo if
    
    
    
    
    # Cambiar el estado basado en los dias que paso en ese mismo
    for (i in 1:length(V(Net)$Days)) {
      if(length(StateConfig) == 3) { # Si es un modelo SIR
        if(V(Net)$Days[i] == S+1) {V(Net)$State[i] <- "B"} else if (V(Net)$Days[i] == I+1) {V(Net)$State[i] <- "C"}
        
        
      } else if(length(StateConfig) == 4) { # Si es un modelo SEIR
        if(V(Net)$Days[i] == S+1) {V(Net)$State[i] <- "B"} else if (V(Net)$Days[i] == E+1) {V(Net)$State[i] <- "C"} else if (V(Net)$Days[i] == I+1) {V(Net)$State[i] <- "D"}
      }
    }
  
    
    
    # Añadir un dia a todos aquellos que iniciaron el proceso infeccioso
    for (i in 1:length(V(Net)$Days)) {
      if(V(Net)$Days[i] > 0) {V(Net)$Days[i] <- V(Net)$Days[i] + 1}
    }
    
    
    
    ####################
    # HAZ UN PLOT AQUI #
    ####################
    my_color <- Colors[as.numeric(as.factor(V(Net)$State))]
    
    plot(Net, layout = LO2, vertex.color = my_color, vertex.size = degree(Adj, V(Adj), "in")*2, 
         vertex.label = "", edge.arrow.size = 0.2, edge.size = 19, 
         vertex.frame.color="white", edge.color = "#E5E3C9")
    
    legend(x=-1.5, y=-1.1, paste(Lol), pch=21, 
           col="lightgray", pt.bg = Colors, pt.cex=2, cex=1.5, bty="n", ncol=1, title = paste(PlotName, "Day", d))
  
  
  #######################################
  # Actualizar la base de datos por dia #
  #######################################
  
  if (length(StateConfig) == 3){
    new_row <- c(as.numeric(table(V(Net)$State == "A")["TRUE"]),
                 as.numeric(table(V(Net)$State == "B")["TRUE"]),
                 as.numeric(table(V(Net)$State == "C")["TRUE"]),
                 d)
    
    Summary <- rbind(Summary, new_row)
  } else if (length(StateConfig) == 4) {
    new_row <- c(as.numeric(table(V(Net)$State == "A")["TRUE"]),
                 as.numeric(table(V(Net)$State == "B")["TRUE"]),
                 as.numeric(table(V(Net)$State == "C")["TRUE"]),
                 as.numeric(table(V(Net)$State == "D")["TRUE"]),
                 d)
    
    Summary <- rbind(Summary, new_row)
  }
  
  d <- d+1 
  
}
return(Summary)
}

Modelo SIR

#### SIR

# Plots de uno por uno
a <- Infectious.Network(Adj, c("S", "I", "R"), c(0, 5, 0), 24, 2, c("#9EA1D4", "#A8D1D1", "#F1F7B5"), 40)

Gift de SIR

# Para hacer un gift
saveGIF(Infectious.Network(Adj, c("S", "I", "R"), c(0, 5, 0), 24, 2, c("#9EA1D4", "#A8D1D1", "#F1F7B5"), 40) ,
        # Nombre del gif
        movie.name = "SIR.gif",
        # Dimensiones
        ani.width  = 1600,
        ani.height = 900,
        # Tiempo de duración de cada frame (segundos)
        interval = 0.2
)
## Output at: SIR.gif
## [1] TRUE

Insertar el gift previamente creado

if (knitr:::is_latex_output()) {
  knitr::asis_output('\\url{https://lh3.googleusercontent.com/pw/AMWts8BbENUtIQFG2YbCKmGyYnrr7siJtiG7QnQYLB5ZSb6PWh0BVIZBOIt4COMqhNSC-faweC1riv3XxS9cxa_y50md7aihQ7PVwM1LQwbQ0PT5nGnKET0TzC8w4huz6I-zRZH0vtTU9-4VXPBGfHEcvpGH=s350-no?authuser=0}')
} else {
  knitr::include_graphics("SIR.gif")
}

Modelo SEIR

# Plots de uno por uno
b <- Infectious.Network(Adj, c("S", "E", "I", "R"), c(0, 2, 5, 0), 24, 2, c("#9EA1D4", "#FD8A8A", "#A8D1D1", "#F1F7B5"), 60)

Gift de SEIR

# Para hacer un gift
saveGIF(Infectious.Network(Adj, c("S", "E", "I", "R"), c(0, 2, 5, 0), 24, 2, c("#9EA1D4", "#FD8A8A", "#A8D1D1", "#F1F7B5"), 60),
        # Nombre del gift
        movie.name = "SEIR.gif",
        # Dimensiones
        ani.width  = 1600,
        ani.height = 900,
        # Tiempo de duración de cada imagen (segundos)
        interval = 0.2
)
## Output at: SEIR.gif
## [1] TRUE

Insertar el gift previamente creado

if (knitr:::is_latex_output()) {
  knitr::asis_output('\\url{https://lh3.googleusercontent.com/pw/AMWts8BbENUtIQFG2YbCKmGyYnrr7siJtiG7QnQYLB5ZSb6PWh0BVIZBOIt4COMqhNSC-faweC1riv3XxS9cxa_y50md7aihQ7PVwM1LQwbQ0PT5nGnKET0TzC8w4huz6I-zRZH0vtTU9-4VXPBGfHEcvpGH=s350-no?authuser=0}')
} else {
  knitr::include_graphics("SEIR.gif")
}

Modelos por compartimentos

Modelo SIR

SIR <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), { 
    
    dS <- -beta*S*I/(S+I+R) 
    dI <- beta*S*I/(S+I+R) -gama*I
    dR <- gama*I 
    
    list(c(dS, dI, dR))
  })
}

###  Condiciones iniciales:

pars <- c(beta = 6, gama = 3)  
condiciones_iniciales <- c(S = 499, I = 1, R = 0)
tiempo <- seq(0, 10, by = 0.0001)
out1 <- ode(condiciones_iniciales, tiempo, SIR, pars)


### Grafica:

matplot(out1[ , 1], out1[ , 2:4], type = "l", xlab = "Time", ylab = "Population", main = "SIR", lwd = 3, col = c("#9EA1D4", "#A8D1D1", "#F1F7B5"))
legend("topright", c("Suceptible", "Infectado / Infected","Recuperado / Recovered"), col = c("#9EA1D4", "#A8D1D1", "#F1F7B5"), lty = 1:4, cex = 0.5)

Modelo SEIR

 SEIR <- function (t, state, parameters) {  
  with (as.list (c (state, parameters)), { 
    dS <- - beta * S * I / (S + E + I + R)  
    dE <- beta * S * I / (S + E + I + R) - delta * E   
    dI <- delta * E - gama * I  
    dR <- gama * I   
    list (c (dS, dE, dI, dR))  
  })
} 

pars <- c (beta = 10, delta = 4, gama = 1) 
condiciones_iniciales <- c (S = 497, E = 2, I = 1, R = 0)  
tiempo <- seq (0, 10, by = 0.001)  
out2 <- ode (condiciones_iniciales, tiempo, SEIR, pars) 

matplot (out2[ , 1], out2[ , 2 : 5], type = "l", xlab = "Time", ylab = "Population", main = "SEIR", lwd = 3, col = c("#9EA1D4", "#FD8A8A", "#A8D1D1", "#F1F7B5")) 
legend ("topright", c ("Susceptible", "Expuesto / Exposed", "Infectado / Infected", "Recuperado / Recovered"), col = c("#9EA1D4", "#FD8A8A", "#A8D1D1", "#F1F7B5"), lty = 1 : 4, cex = 0.5)

Comparación de la función de Redes y desolve

Quitar NA’s de los resultados de las funciones (NA=0)

# Tablas del resultado de la funcion de redes

# Quitar las NA's de las tablas (NA=0)

a[is.na(a)] <- 0


b[is.na(b)] <- 0

Modelo SIR

Función de redes

# Funcion de redes 

matplot(a[ -1, 4], a[ -1, 1:3], type = "l", xlab = "Time", ylab = "Population", main = "SIR con función de redes", lwd = 3, col = c("#9EA1D4", "#A8D1D1", "#F1F7B5"))
legend("topright", c("Suceptible", "Infectado / Infected","Recuperado / Recovered"), col = c("#9EA1D4", "#A8D1D1", "#F1F7B5"), lty = 1:4, cex = 0.5)

Desolve

# Desolve

matplot(out1[ , 1], out1[ , 2:4], type = "l", xlab = "Time", ylab = "Population", main = "SIR con desolve", lwd = 3, col = c("#9EA1D4", "#A8D1D1", "#F1F7B5"))
legend("topright", c("Suceptible", "Infectado / Infected","Recuperado / Recovered"), col = c("#9EA1D4", "#A8D1D1", "#F1F7B5"), lty = 1:4, cex = 0.5)

Modelo SEIR

# Funcion de redes

matplot (b[ -1, 5], b[ -1, 1:4], type = "l", xlab = "Time", ylab = "Population", main = "SEIR con función de redes", lwd = 3, col = c("#9EA1D4", "#FD8A8A", "#A8D1D1", "#F1F7B5")) 
legend ("topright", c ("Susceptible", "Expuesto / Exposed", "Infectado / Infected", "Recuperado / Recovered"), col = c("#9EA1D4", "#FD8A8A", "#A8D1D1", "#F1F7B5"), lty = 1 : 4, cex = 0.5) 

# Desolve 

matplot (out2[ , 1], out2[ , 2 : 5], type = "l", xlab = "Time", ylab = "Population", main = "SEIR con desolve", lwd = 3, col = c("#9EA1D4", "#FD8A8A", "#A8D1D1", "#F1F7B5")) 
legend ("topright", c ("Susceptible", "Expuesto / Exposed", "Infectado / Infected", "Recuperado / Recovered"), col = c("#9EA1D4", "#FD8A8A", "#A8D1D1", "#F1F7B5"), lty = 1 : 4, cex = 0.5)

Conclusiones y perspectivas

La función de redes, mostró resultados óptimos comparado con la función ode de desolve, ambos métodos llegaron a un punto de equilibrio.

Comparaciones entre los métodos:

En el modelo SIR: El pico de personas infectadas en la función de redes fue mayor que en desolve, los susceptibles se redujeron en menos tiempo en desolve y permanecieron más personas susceptibles. No se observa una grande entre los recuperados.

En el modelo SEIR: En la función de redes el pico de infectados y expuestos fue menor que en desolve y en un mayor tiempo. La cantidad de personas recuperadas es mayor en desolve.

La mayor diferencia entre ambos métodos se observó en el tiempo. Sin embargo, en desolve el tiempo = 10; en este caso, es relativo, ya que no se especifica que sean 10 días y pueden ser 10 semanas, etc. En cambio, en la función de redes sí se basa en la cantidad de días solamente. Como perspectivas, se puede buscar implementar la función a modelos con demografía y a otros modelos con otros estados, como crónico, en tratamiento, etc.

Referencias